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pair of polymers to fully overlap of order 2/c^T, where 
ks is Boltzmann's constant and T is the temperature 
[T8H28] . Dendrimers, due to the nature of their chem- 
ical architecture, can have an effective interaction with 
a higher value of n; such systems form so-called 'cluster 
crystals' [28 and there has been a great deal of interest 
in soft potential models for these systems [29H43] . 

Here we consider a model two-dimensional system of 
soft 'particles' that interact via the potential 
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Formation of quasicrystals and other exotic crystal structures in a soft-core fluid: 
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A two-dimensional system of soft particles interacting via a two- length-scale potential is studied. 
Density functional theory reveals the existence of a fluid phase and two crystalline phases with 
different lattice spacing. Of these the larger lattice spacing crystalline phase takes the form of a 
'crystal liquid' state with a small fraction of mobile particles. Near the transition between this phase 
and the smaller lattice spacing crystalline phase, quasicrystalline structures may be created by a 
competition between linear instability at one scale and nonlinear selection of the other. The results 
are confirmed using Brownian dynamics simulations. 
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Regular crystals are ordered arrangements of atoms 
or molecules with rotation and translation symmetries. 
Quasicrystals, discovered in 1984 pQ, lack the lattice sym- 
metries of crystals and yet have discrete Fourier spectra. 
Quasicrystals have been found not only in metals but 
also in systems of nanoparticles [2 , mesoporous silica [3] , 
and soft-matter systems [4]. The latter include micellar 
melts [5j [6] formed, e.g., from linear, dendrimer or star 
block copolymers, comprising a stiff hydrophobic poly- 
mer core surrounded by a corona of flexible hydrophilic 
polymer chains. The main theoretical approach to inves- 
tigating the stability of metallic or micellar quasicrystals 
involves minimising an appropriate energy, but the prin- 
ciple underlying their stability is not known [7J [8] . 

Patterns with quasicrystalline structure, or quasipat- 
terns, were discovered in Faraday wave experiments in 
the 1990s [9j[10]. Many quasipatterns are characterised 
by two length scales [TTJ [12] , and recent work suggests 
that nonlinear interaction between these scales can sta- 
bilise such patterns p~3j [14] and this is so for soft-matter 
quasicrystals as well [T5HT7] . In this Letter, we explore 
this mechanism in detail, by considering a system of par- 
ticles interacting via a simple potential with two length 
scales. We compute the phase diagram, finding crystals 
of each length scale and a region of metastable quasicrys- 
tals with both length scales prominent in their Fourier 
spectra. These form in a region where the linear growth 
of density fluctuations in a quenched uniform fluid favors 
one length scale but nonlinear stability favors the other. 

The effective 'coarse-grained' interaction potentials be- 
tween the centres of mass of polymers, dendrimers or 
other such macromolecules, are soft. By this we mean 
that they are finite for all separation distances r, be- 
cause the centre of mass of such soft objects does not 
necessarily coincide with any individual monomer. The 
soft effective pair potential between such particles can be 
approximated as V(r) = ee~( r / R ^ n . Simple linear poly- 
mers correspond to the case n = 2 with the length R of 
order the radius of gyration and the energy e for such a 



(i) 



The energy for complete overlap is (1 + A)e. The most 
important feature of this potential is that it has a 'shoul- 
der' when the dimensionless parameter i/0, with two 
length scales. The radius of the core is R and the ra- 
dius of the shoulder is R s > R. Such a potential is a 
simple coarse-grained model for the effective interaction 
between polymers or micelles formed, e.g., from linear, 
dendrimer or star block copolymers, which have a stiff 
hydrophobic core surrounded by a corona of flexible hy- 
drophilic chains. A related, piecewise constant poten- 
tial is used in Ref. [17] • In the following we set the di- 
mensionless interaction energy parameter f3e — 1, where 
j3 = (&bT) _1 , and fix the ratio of the two length scales 
to be R s /R = 1.855. 

We use density functional theory (DFT) f44U46] to 
study this system. The grand free energy is 

Q[p(r)] = k B T J drp(r)[ln A 2 p(r)-l] 

+ 7- ex [p(r)] + ydr($(r)-/x)p(r), (2) 

which is a functional of the one-body (number) density of 
the particles, p(r), where r = (x, y). The first term is the 
ideal-gas contribution to the free energy, A is the (irrel- 
evant) thermal de Broglie wavelength, \i is the chemical 
potential, <£>(r) is any external potential that may be con- 
fining the system and F ex [p(r)] is the excess Helmholtz 



free energy from the interactions between the particles. 
The equilibrium density distribution is that which min- 
imises Q[p(r)]; the corresponding minimum is the ther- 
modynamic grand potential of the system. For a system 
in the bulk fluid state (i.e., where 3>(r) = 0), the min- 
imising density is uniform, p = po. However, for other 
state points, when the system freezes to form a solid, Q is 
minimised by nonuniform density distributions, exhibit- 
ing sharp peaks. For the systems of soft-core particles 
considered here, one may approximate T ex as [18] : 



?exW)\ = lJdrJdr'p(r)V(\r-r'\)p(r'). 



(3) 



This functional generates the random phase approxi- 
mation (RPA) for the pair direct correlation function 
c (2) (r y) ^ -p **f* (rl) = -f3V(\r-v>\) IMSS]. If we 



6p(r)6p(r') 

assume that these are Brownian particles with dynamics 
governed by 



-TV i U({r i },t) + TX i (t), 



(4) 



where the index i = 1,...,7V labels the particles, 

U({ri},t) = £^1 $( r + £i# V ( r i ~ t j) is the P° ten " 
tial energy of the system and X^ (t) is a white noise term, 
we can investigate the dynamics of the system using Dy- 
namic Density Functional Theory (DDFT) [47H50] in the 
form 



dp(r,t) 
dt 



rv 



p(r,t)V 



Sn[p(r,t)\ 



6p(r,t) 



(5) 



where p(r,i) is now the time-dependent nonequilibrium 
one-body density profile and T = (3D is the mobility. 
Here D is the diffusion coefficient. In deriving this the- 
ory we have used the equilibrium Helmholtz free energy 
functional T to approximate the unknown nonequilib- 
rium free energy. 

Fig. [IJa) shows the equilibrium phase diagram calcu- 
lated using Picard iteration of the DFT Euler-Lagrange 
equation, starting either from the profile for a nearby 
state point or a uniform density profile with a small 
random value added to each point. As the fluid density 
is increased, the system freezes to form one of two 
distinct solid phases (Fig. [2]): for larger values of A the 
system forms a hexagonal crystal with a large lattice 
spacing (referred to as 'crystal A'), but for smaller 
values of A it forms a hexagonal crystal with a much 
smaller lattice spacing ('crystal B'). The red regions in 
the diagram denote thermodynamic coexistence between 
two different phases at the same temperature, pressure 
and chemical potential. 

To understand the phase diagram we study the struc- 
ture and stability of a uniform liquid with density po 
and $(r) = 0. We follow 031 S3 ED [52] and expand 
Eq. ([5]) in powers of p(r, t) = p(r, t) — po. Retain- 
ing only linear terms, we find that the growth/decay 
of different Fourier modes of wave number k follows 



p(k,t) = p(fc, 0)exp[u;(fc)£], where uo(k) satisfies the dis- 
persion relation [49, 52 



j(k) = -Tk B Tk 2 (l-p c(k)). 



(6) 



Here c(k) is the Fourier transform of the pair direct cor- 
relation function; within RPA c(k) = —(3V(k), where 
V (fc) is the Fourier transform of the pair potential in Eq. 
([l]). In an equilibrium fluid the static structure factor 
S(k) = (1 — pod(k))~ 1 > for all values of k; such a fluid 
is therefore stable [53] . Within RPA the two length scales 
in the pair potential lead, for certain ranges of parameter 
values, to a static structure factor S(k) with two peaks. 
Fig. mb)-(d) shows that as A increases the smaller wave 
number peak in S(k) grows and comes to dominate the 
larger wave number peak. Fig. [IJe)-(g) shows analo- 
gous behavior of the dispersion relation at several points 
in or on the boundary of the linearly unstable region 
u(k max ) = where k max is the wave number of the higher 
peak (blue dashed line in Fig. []Ja)): as A increases the 
instability shifts from large wave number (Fig. [Tig)) to 
small wave number (Fig. file)). The short and long length 
scales are simultaneously marginally stable at A = 1.067 
and poR 2 = 2.95 (Fig. [Iff)); this point lies on the pink 
dotted line in Fig. [TJa) corresponding to a pair of equal 
height peaks in the dispersion relation. Above (below) 
this line, the peak at smaller (larger) wave number k is 
higher, indicating that the longer (shorter) length scale 
density fluctuations grow the fastest. The black double 
dotted lines indicate the location of oj(k max ) = for the 
lower peak in the dispersion relation. When the system is 
quenched from a stable liquid state to a state point with 
density po above the blue dashed line, certain wave num- 
bers will grow as described by the dispersion relation (|6| . 
Fig. [2] third panel in the top row, shows the density 
profile of the larger lattice spacing crystal A phase for a 
state point not far from the transition to the smaller lat- 
tice spacing crystal B phase. However, the panel below 
displaying \og[R 2 p(r)] reveals an interconnected network 
of 'channels' between the density peaks. The particles 
contributing to this part of the density profile are fluid 
in the sense that they can move freely throughout the 
whole system, unlike the majority of the particles that 
are located in the density peaks fixed on the lattice 
sites. We refer to this structure as the 'crystal liquid' 
state. This state minimises the free energy for A > A COl 
where A co is the value at coexistence, and is reminiscent 
of the way electrons arrange themselves in a metal, 
where the vast majority are bound on the lattice sites 
to the metal atomic nuclei. However, there are some, 
the conduction band electrons, that are free to move 
throughout the crystal. The crucial differences between 
the present system and a metal are that (i) the present 
particles are classical, and (ii) they self-organise to form 
this remarkable structure - there is no equivalent of the 
atomic nuclei to provide a periodic array of potential 
wells into which the majority of the particles condense. 
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FIG. 1: (color online) Phase diagram, static structure factor S(k) and dispersion relation cj(k) for f3e — 1 and R s /R — 1.855. 
(a) The bulk system phase diagram in the (poR 2 , A) plane. The system exhibits a uniform fluid phase and two crystal phases: 
the larger lattice spacing 'crystal A' phase and the smaller lattice spacing 'crystal B' phase. The regions filled in red denote 
areas where there is two-phase coexistence between the different phases. The blue dashed line denotes the linear instability 
threshold for the liquid phase while the pink dotted line terminating in a circle is the locus where the two peaks in the 
dispersion relation (|6| have the same height. The circle denotes the point where the smaller-/c peak disappears, (b)-(d) S(k) 
for (b) (p R 2 , A) = (0.8, 1.5), (c) (1.8, 1.067), (d) (2, 0.7).(e)-(g) u(k) for (e) (p R 2 , A) = (2, 1.5), (f) (2.95, 1.067), (g) (3.5, 0.7). 
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FIG. 2: (color online) Density profiles from DFT showing R 2 p(r) (upper panels) and \og[R 2 p(r)] (lower panels). From left to 
right: (poR 2 , A) = (4.4, 0.7) (typical of the small length scale crystal 'crystal B'), (3.5, 0.76), (4.0, 0.8) (both near the transition 
from 'crystal A' to 'crystal B') and (2.7, 2) (typical of the large length scale crystal 'crystal A'). The second set of panels shows 
quasicrystalline ordering with numerous defects, while the third set reveals a network of connected density, indicating that 
the particles in this part of the crystal are fluid, and able to move throughout the system. There are also similar connected 
fragments in the disordered (3.5,0.76) profile, but because of the disorder, these do not percolate the system. 



As A decreases towards A co1 the fraction of particles 
in this 'conduction band' increases and in fact it is 
this growing number that triggers the formation of the 
smaller length scale crystal: these mobile particles freeze 
to form the extra peaks of crystal B. 



To confirm the existence of the 'crystal liquid' state 
we calculated the density profile for a system confined 
within a square confining potential <£> of size L x L with 
hard walls, and compared the results with Brownian dy- 
namics (BD) simulations, i.e., simulations of N particles 
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FIG. 3: log[i? 2 p(r)] for a system of N = 600 particles with 
(poR 2 , A) = (4.0, 0.8) confined in a square region of side L = 
10R. The left panel displays the BD simulation result, while 
the right panel is from DFT. The system forms crystal A with 
a density profile consisting of an array of peaks surrounded 
by a connected network within which the particles are free to 
move - the 'crystal liquid' phase. 
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FIG. 5: Grand potential Q, per unit area for f5fi — 39 as a func- 
tion of A for the two different crystal structures together with 
the free energy for the quasicrystalline solution displayed in 
Fig. [4] There is a point where all three have almost the same 
free energy, but the quasicrystal solution never corresponds to 
the global minimum of the free energy. The crystal A phase 
is of 'crystal liquid' type throughout the range of A shown. 
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FIG. 4: Left: \og[R 2 p(r)} from DFT, for (p R 2 ,A) = 
(3.5,0.8). Right: the corresponding Fourier transforms. The 
12-fold symmetry is indicative of quasicrystalline ordering. 
The upper density profile was obtained from random initial 
conditions, while the lower one was started from initial con- 
ditions with quasicrystalline symmetry. 



evolving according to Eq. Q. Averaging over the posi- 
tions of the particles to calculate the density profile, we 
find remarkably good agreement between the DFT and 
the BD results (Fig. [3]). 

Formation of a c quasicrystaV: A striking aspect of the 
phase diagram in Fig. HJ a) is that the phase transition 
between the two different crystal phases (thin red region) 
is well away from where the two peaks in the dispersion 
relation have the same height (pink dotted line). A uni- 
form system quenched to the region above the coexistence 
of the two crystal phases but below this line will initially 
generate small length scale density fluctuations and the 



system behaves as if it were going to form crystal B. How- 
ever, the true minimum of the free energy is the larger 
length scale crystal. Thus, as growing density fluctua- 
tions reach the nonlinear regime, the system seeks to go 
to the longer length scale structure but the smaller length 
scale imprinted from the linear growth regime leads to 
frustration. Sometimes the system is able to evolve to the 
larger length scale crystal; at other times it stays stuck 
in the metastable small length scale crystal B structure. 
However, often the system forms a state with density 
peaks on both length scales, but no long range order. 
In Fig. [4] we display two rather striking density profiles 
calculated at a state point in this region of the phase di- 
agram. The upper density profile was calculated using 
Picard iteration starting from random initial conditions. 
The density profile has many defects, but it has definite 
quasicrystalline ordering, as can be seen from the corre- 
sponding Fourier transform. The lower panels in Fig. [4] 
show a defect-free quasicrystal approximant, started from 
carefully chosen initial conditions. The two wavenumbers 
k\R = 3.2 and fei? = 6.0 corresponding to the maxima 
in uj(k) are indicated in the Fourier transforms. 

The Picard iteration of the Euler-Lagrange equation 
corresponds to fictitious dynamics since it does not con- 
serve the total number of particles in the system, N = 
Jdrp(r). The true dynamics is governed by the DDFT 
equation ^. Evolving this equation can be rather slow 
but in most cases the same qualitative behavior is ob- 
served, although the system sometimes gets stuck in the 
smaller length scale crystal B that initially forms in the 
linear growth regime. This is a consequence of the con- 
served dynamics. For /3e = 1, R s /R = 1.855 the qua- 
sicrystals we find are never the minimum free energy state 
(Fig. [5]). The quasicrystaline state in Fig. [4] remains sta- 
ble for 1.77 < R s /R < 2.18, but we have not calculated 
the full phase diagram for R s ^ 1.855. We believe it may 
be possible to use nonlinear dynamics techniques [14] to 
compute the stability properties of these states by re- 
ducing the DDFT description in Eq. (pi to a phase field 



crystal model, cf. [52j [54H57] . We expect that the ob- 
served quasicrystal formation mechanism (linear growth 
of one length scale, but nonlinear selection favoring an- 
other) may well apply more generally. 
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